rm(list=ls())
setwd("folder/Analysis")
library(haven)
library(dplyr)
library(tidyverse)
library(haven)
library(foreign)
library(ggplot2)

all_data <- read.csv("taiwandata.csv")

## set reform year
all_data <- all_data %>%
  mutate(PostReform = case_when(
    Year == "2006" ~ "0", 
    Year == "2010" ~ "0",
    Year =="2014"  ~ "0",
    Year =="2018"  ~ "0",
    Year =="2021"  ~ "1"
  ))

all_data <- na.omit(all_data)

all_data$TrustCourt <- as.numeric(all_data$TrustCourt)
all_data$AgeGroup <- as.numeric(all_data$AgeGroup)
all_data$TrustExe <- as.numeric(all_data$TrustExe)
all_data$TrustLegis <- as.numeric(all_data$TrustLegis)
all_data$engagement <-as.factor(all_data$engagement)

####
fmla <- TrustCourt ~ PostReform + MediaFreq + PostReform*MediaFreq + 
  Male +  AgeGroup + college + LikeKMT  + TrustLegis + engagement

fit <- lm(fmla,data=all_data)

## matching
library(MatchIt)
library(sandwich)
library(lmtest)

all_data$treated <- ifelse(all_data$MediaFreq == 1, "1", "0")


match1 <- matchit(treated ~  AgeGroup + college + LikeKMT,
                  method = "nearest", 
                  distance = "mahalanobis", 
                  data=all_data)
summary(match1, covariates = T)

match1.data <- match.data(match1)


Original <- lm(fmla,data=all_data)
summary(Original)

Matched <- lm(fmla,data=match1.data, weights = weights)
summary(Matched)

sd_dv <- sd(match1.data$TrustCourt)

##
install.packages("stargazer")
library(stargazer)

stargazer(Original, Matched)

##
##
## Placebo Test
##
##  create a fake dataset wtih data from pre-reform
fake2006_match1.data <- match1.data %>% 
  mutate(PostReform = case_when(
    Year == "2006" ~ "1",
    Year == "2010" ~ "0",
    Year =="2014"  ~ "0",
    Year =="2018"  ~ "0"
  ))


fake2010_match1.data <- match1.data %>% 
  mutate(PostReform = case_when(
    Year == "2006" ~ "0",
    Year == "2010" ~ "1",
    Year =="2014"  ~ "0",
    Year =="2018"  ~ "0"
  ))

fake2014_match1.data <- match1.data %>% 
  mutate(PostReform = case_when(
    Year == "2006" ~ "0",
    Year == "2010" ~ "0",
    Year =="2014"  ~ "1",
    Year =="2018"  ~ "0"
  ))

fake2018_match1.data <- match1.data %>% 
  mutate(PostReform = case_when(
    Year == "2006" ~ "0",
    Year == "2010" ~ "0",
    Year =="2014"  ~ "0",
    Year =="2018"  ~ "1"
  ))


fake2006 <- na.omit(fake2006_match1.data)
fake2010 <- na.omit(fake2010_match1.data)
fake2014 <- na.omit(fake2014_match1.data)
fake2018 <- na.omit(fake2018_match1.data)

##

fake2006$TrustCourt <- as.numeric(fake2006$TrustCourt)
fake2006$AgeGroup <- as.numeric(fake2006$AgeGroup)
fake2006$TrustLegis <- as.numeric(fake2006$TrustLegis)


fake2010$TrustCourt <- as.numeric(fake2010$TrustCourt)
fake2010$AgeGroup <- as.numeric(fake2010$AgeGroup)
fake2010$TrustLegis <- as.numeric(fake2010$TrustLegis)
####
fake2014$TrustCourt <- as.numeric(fake2014$TrustCourt)
fake2014$AgeGroup <- as.numeric(fake2014$AgeGroup)
fake2014$TrustLegis <- as.numeric(fake2014$TrustLegis)
####
fake2018$TrustCourt <- as.numeric(fake2018$TrustCourt)
fake2018$AgeGroup <- as.numeric(fake2018$AgeGroup)
fake2018$TrustLegis <- as.numeric(fake2018$TrustLegis)

##
real <- lm(fmla,data=match1.data)

fake1 <- lm(fmla,data=fake2006)
fake2 <- lm(fmla,data=fake2010)
fake3 <- lm(fmla,data=fake2014)
fake4 <- lm(fmla,data=fake2018)

##
summary(fake1)
summary(fake2)
summary(fake3)
summary(fake4)
summary(real)

#####
library(MCMCpack)
install.packages("matrixStats")
library(matrixStats)
library(Matrix)
## Simulate coefficients
beta2021 <- mvrnorm(n=1000, mu=coef(real), Sigma=vcov(real))

####
beta2006 <- mvrnorm(n=1000, mu=coef(fake1), Sigma=vcov(fake1))
beta2010 <- mvrnorm(n=1000, mu=coef(fake2), Sigma=vcov(fake2))
beta2014 <- mvrnorm(n=1000, mu=coef(fake3), Sigma=vcov(fake3))
beta2018 <- mvrnorm(n=1000, mu=coef(fake4), Sigma=vcov(fake4))


effect2006 <-beta2006[,11]
effect2006 <- as.data.frame(effect2006)

effect2010 <-beta2010[,11]
effect2010 <- as.data.frame(effect2010)

###
effect2014 <-beta2014[,11]
effect2014 <- as.data.frame(effect2014)

##
effect2018 <-beta2018[,11]
effect2018 <- as.data.frame(effect2018)

effect2021 <- beta2021[,11]
effect2021 <- as.data.frame(effect2021)

fakes <- cbind(effect2006, effect2010, effect2014, effect2018,effect2021)
fakes<-as.matrix(fakes)

fakes_mean <- colMeans(fakes)
fakes_q975 <- colQuantiles(fakes, probs=.975)
fakes_q025 <- colQuantiles(fakes, probs=.025)

fakes_cbind <- cbind(fakes_q025, fakes_mean, fakes_q975)

colnames(fakes_cbind) <- c("Lower", "Mean", "Upper")
fakes_cbind <- as.data.frame(fakes_cbind)
fakes_cbind$Year <- c("2006","2010", "2014", "2018", "2021")
fakes_cbind$Year <-as.numeric(fakes_cbind$Year)
### 
fig1 <- ggplot(fakes_cbind, aes(x = Year, y = Mean)) +
  geom_pointrange(aes(ymin = Lower, ymax = Upper), size = 1, 
                  position = position_dodge(width = 0.5))  +
  ylim(c(-0.2, 0.3))+
  geom_hline(yintercept=0, linetype="dashed", color = "black") +
  geom_vline(xintercept = 2019, linetype="dashed", color = "red", size = 1) +
  scale_x_continuous(breaks = c(2006:2021), 
                     labels = c(2006:2021)) +
  labs(title="", size = 0.5) +
  ylab("Estimated Coefficients On The Treatment Term") + xlab("Hypothetical Adoption Year of Judicial Reform") + 
  theme_bw() + theme(legend.position="bottom") 
fig1
